# -*- coding: utf-8 -*-
"""
Created on Fri Oct 25 22:19:23 2024
设计了二维矩形房间内的势流流动模型
此模型关注房间内整体的闭合流线及流动特性，包括源、汇，属于弯折拐角壁面流动bended-corner wall flow的研究

@author: wangduo
"""
import numpy as np
import matplotlib.pyplot as plt
import math
from geometry import Grid
from visualize import CornerFlowPlotter, SosiFlowPlotter
import scipy.special as sp
from scipy.optimize import root_scalar
from numpy.polynomial import Polynomial
# embed figures in the notebook
#%matplotlib inline
#%config InlineBackend.figure_format='svg'# 让图片看起来像矢量图，更清晰

# 椭圆函数求解器类：根据纵横比计算出SC变换的参数k
class ParameterSolver:
    def __init__(self, L, H):
        """
        初始化参数求解器，输入 L 和 H
        L: 长度参数
        H: 高度参数
        """
        self.L = L
        self.H = H
        self.k_value = None  # 模数 k
        self.k_prime_value = None  # 补模数 k'
        self.C_value = None  # 参数 C
    
    def solve_k_and_k_prime(self):
        """
        通过方程 (1/2 * L / H) * K(k') - K(k) = 0 求解模数 k 和补模数 k'
        """
        # 定义目标方程
        def equation(k):
            k_prime = np.sqrt(1 - k**2)  # 补模数 k'
            K_k = sp.ellipk(k**2)        # 计算 K(k)
            K_k_prime = sp.ellipk(k_prime**2)  # 计算 K(k')
            
            # 方程 (1/2 * L / H) * K(k') - K(k) = 0
            return (1/2 * self.L / self.H) * K_k_prime - K_k
        
        # 使用数值方法求解 k
        sol = root_scalar(equation, bracket=[0, 1], method='bisect')
        
        if sol.converged:
            self.k_value = sol.root
            self.k_prime_value = np.sqrt(1 - self.k_value**2)
        else:
            raise RuntimeError("求解模数 k 和补模数 k' 时未收敛")

    def verify_k_and_k_prime(self):
        """
        验证 k^2 + k'^2 是否等于 1
        """
        if self.k_value is None or self.k_prime_value is None:
            raise ValueError("请先求解 k 和 k'")
        
        sum_of_squares = self.k_value**2 + self.k_prime_value**2
        if not np.isclose(sum_of_squares, 1):
            raise ValueError(f"验证失败：k^2 + k'^2 = {sum_of_squares} 不等于 1")
        else:
            print("验证成功：k^2 + k'^2 接近 1")

    def compute_C(self):
        """
        通过方程 1/C - K(k') = 0 计算 C
        """
        if self.k_prime_value is None:
            raise ValueError("请先求解 k 和 k'")
        
        K_k_prime = sp.ellipk(self.k_prime_value**2)  # 计算 K(k')
        self.C_value = 1 / K_k_prime
        return self.C_value
    
    def solve(self):
        """
        执行求解过程，输出 C、k、k' 并进行验证
        """
        print("开始求解模数 k 和补模数 k'...")
        self.solve_k_and_k_prime()  # 求解 k 和 k'
        
        print(f"求解得到的模数 k: {self.k_value}")
        print(f"求解得到的补模数 k': {self.k_prime_value}")
        
        # 验证 k^2 + k'^2 是否等于 1
        print("验证 k 和 k'...")
        self.verify_k_and_k_prime()
        
        # 计算 C
        print("计算参数 C...")
        C = self.compute_C()
        print(f"通过方程 1/C - K(k') = 0 计算得到的 C: {C}")
        
        return self.C_value, self.k_value, self.k_prime_value
    

# 室内源汇类
class RoomSourceSink:
    """
    Contains information related to a source and a sink in the room.
    """
    def __init__(self, length, height, length_1, height_1, length_2, height_2, strength, C, k):
        """
        Sets the location and strength of the singularity.
        
        Parameters
        ----------
        length: float
            length of room (m)
        height: float
            height of room (m)
        length_1: float
            distance from the source to the left wall (m)
        height_1: float
            height of the source from the floor (m)
        length_2: float
            distance from the sink to the left wall (m)
        height_2: float
            height of the sink from the floor (m)
        strength: float
            Strength of source and sink.
        C: float
            parameter for the SC transformation function
        k: float
            modulus number for elliptic functions
        """
        self.Q = strength  # 源、汇强度
        self.L = length
        self.H = height
        self.L1 = length_1
        self.H1 = height_1
        self.L2 = length_2
        self.H2 = height_2

        # 计算无量纲的源点和汇点坐标，稍后会调用 points_trans
        self.z_so = None  # 源点复数坐标
        self.z_si = None  # 汇点复数坐标
        
        self.C0 = C  # SC 变换的常数
        self.k0 = k  # 模数
        self.m0 = k**2  # 模数平方
        self.k0p = np.sqrt(1 - k**2)  # 补模数
        self.m0p = 1 - self.m0  # 补模数平方
        self.epsilon = 1e-10  # 小数值避免除零

    def compute_z_coordinates(self):
        '''
        计算源点和汇点的无量纲坐标，并保存到类变量 z_so 和 z_si 中
        '''
        # 检查输入合法性
        if self.H <= 0:
            raise ValueError("房间高度 H 必须大于 0")
        if self.L1 < 0 or self.L1 > self.L or self.L2 < 0 or self.L2 > self.L:
            raise ValueError("源点或汇点的水平位置超出了房间宽度范围")
        if self.H1 < 0 or self.H1 > self.H or self.H2 < 0 or self.H2 > self.H:
            raise ValueError("源点或汇点的垂直位置超出了房间高度范围")

        # 有量纲坐标（复数表示）
        E = self.L1 - self.L / 2 + 1j * self.H1    # 源点
        F = self.L2 - self.L / 2 + 1j * self.H2    # 汇点

        # 无量纲坐标
        self.z_so = E / self.H  # 源点无量纲坐标
        self.z_si = F / self.H  # 汇点无量纲坐标

        #print(f"计算得到的源点坐标 z_so: {self.z_so}")
        #print(f"计算得到的汇点坐标 z_si: {self.z_si}")

    # 源汇复势
    def cplx_W(self, X, Y):
        """
        计算复位势及其分量
        :param X: 笛卡尔坐标 X
        :param Y: 笛卡尔坐标 Y
        :return: w_ptn (复势), phi (位势), psi (流函数)
        """
        # 确保先计算源点和汇点的复数坐标
        if self.z_so is None or self.z_si is None:
            self.compute_z_coordinates()
            
        Z = X + Y * 1j  # 复数形式的坐标 Z = X + iY

        # 定义SC变换函数 zeta = sn(z / C0; k**2)
        def SCT_g(X, Y):
            """
            计算复数域的椭圆正弦函数 sn(z, k)
            Parameters:
            X -- 实部
            Y -- 虚部
            Returns:
            复数域的椭圆正弦函数结果 zeta
            """
            # 归一化
            use_X = X / self.C0
            use_Y = Y / self.C0

            # 椭圆函数的计算
            sn_X_k0, cn_X_k0, dn_X_k0, _ = sp.ellipj(use_X, self.m0)
            sn_Y_k0, cn_Y_k0, dn_Y_k0, _ = sp.ellipj(use_Y, self.m0)
            
            sn_X_k0p, cn_X_k0p, dn_X_k0p, _ = sp.ellipj(use_X, self.m0p)
            sn_Y_k0p, cn_Y_k0p, dn_Y_k0p, _ = sp.ellipj(use_Y, self.m0p)

            # 计算分子和分母
            num = sn_X_k0 * dn_Y_k0p + 1j * sn_Y_k0p * cn_Y_k0p * cn_X_k0 * dn_X_k0
            den = cn_Y_k0p**2 + self.k0**2 * sn_X_k0**2 * sn_Y_k0p**2

            zeta = num / den
            return zeta

        # 初始化返回数组
        w_ptn = np.zeros_like(Z, dtype=np.complex128)
        phi = np.zeros_like(Z, dtype=np.float64)
        psi = np.zeros_like(Z, dtype=np.float64)

        zeta1 = SCT_g(np.real(self.z_so), np.imag(self.z_so))  # 源点变换
        zeta2 = SCT_g(np.real(self.z_si), np.imag(self.z_si))  # 汇点变换

        # 遍历 Z 数组
        it = np.nditer(Z, flags=['multi_index'])
        while not it.finished:
            z = it[0]
            zeta = SCT_g(np.real(z), np.imag(z))  # 当前点的 SC 变换

            # 避免分母为零，添加 epsilon
            #denominator_z1 = zeta - zeta1 + self.epsilon
            #denominator_z2 = zeta - zeta2 + self.epsilon

            # 计算复势函数 w_ptn
            w_ptn[it.multi_index] = (self.Q / np.pi) * (np.log(zeta - zeta1) - np.log(zeta - zeta2))

            # 提取位势和流函数
            phi[it.multi_index] = w_ptn[it.multi_index].real
            psi[it.multi_index] = w_ptn[it.multi_index].imag

            it.iternext()  # 移动到下一个坐标点

        return w_ptn, phi, psi

    # 源汇复速度
    def cplx_velocity(self, X, Y):
        """
        计算复势 w 对复坐标 Z 的导数
        :param X: 笛卡尔坐标 X
        :param Y: 笛卡尔坐标 Y
        :return: dw_dZ (复势的导数)
        """
        # 确保先计算源点和汇点的复数坐标
        if self.z_so is None or self.z_si is None:
            self.compute_z_coordinates()
            
        Z = X + Y * 1j  # 复数形式的坐标 Z = X + iY
    
        # 定义SC变换函数 zeta = sn(z / C0; k**2)
        def SCT_g(X, Y):
            """
            计算复数域的椭圆正弦函数 sn(z, k)
            Parameters:
            X -- 实部
            Y -- 虚部
            Returns:
            复数域的椭圆正弦函数结果 zeta
            """
            # 归一化
            use_X = X / self.C0
            use_Y = Y / self.C0

            # 椭圆函数的计算
            sn_X_k0, cn_X_k0, dn_X_k0, _ = sp.ellipj(use_X, self.m0)
            sn_Y_k0, cn_Y_k0, dn_Y_k0, _ = sp.ellipj(use_Y, self.m0)
            
            sn_X_k0p, cn_X_k0p, dn_X_k0p, _ = sp.ellipj(use_X, self.m0p)
            sn_Y_k0p, cn_Y_k0p, dn_Y_k0p, _ = sp.ellipj(use_Y, self.m0p)

            # 计算分子和分母
            num = sn_X_k0 * dn_Y_k0p + 1j * sn_Y_k0p * cn_Y_k0p * cn_X_k0 * dn_X_k0
            den = cn_Y_k0p**2 + self.k0**2 * sn_X_k0**2 * sn_Y_k0p**2

            zeta = num / den
            return zeta
    
        # 椭圆正弦函数的导数 dSCTdZ(z / C0; k**2)
        def SCT_g_prime(X, Y):
            """
            z = x + iy
            计算椭圆正弦函数的导数 cn(z, k) * dn(z, k)
            Z -- 复数形式的坐标
            """
            # step0 归一化
            use_X = X / self.C0
            use_Y = Y / self.C0
            # step1 实部、虚部分离计算模数的椭圆函数
            sn_X_k0, cn_X_k0, dn_X_k0, _ = sp.ellipj(use_X, self.m0)
            sn_Y_k0, cn_Y_k0, dn_Y_k0, _ = sp.ellipj(use_Y, self.m0)
            # step2 实部、虚部分离计算补模数的椭圆函数
            sn_X_k0p, cn_X_k0p, dn_X_k0p, _ = sp.ellipj(use_X, self.m0p)
            sn_Y_k0p, cn_Y_k0p, dn_Y_k0p, _ = sp.ellipj(use_Y, self.m0p)
            # step3 组件计算
            sn_iY_k0 = 1j * (sn_Y_k0p/cn_Y_k0p)
            cn_iY_k0 = 1 / cn_Y_k0p
            dn_iY_k0 = dn_Y_k0p/cn_Y_k0p

            CN_num = sn_X_k0*cn_X_k0*dn_iY_k0-sn_iY_k0*cn_iY_k0*dn_X_k0
            CN_den = sn_X_k0*cn_iY_k0*dn_iY_k0-sn_iY_k0*cn_X_k0*dn_X_k0
            CN_Z_k0 = CN_num / CN_den

            DN_num = dn_X_k0*cn_Y_k0p*dn_Y_k0p-1j*self.m0*sn_X_k0*cn_X_k0*sn_Y_k0p
            DN_den = cn_Y_k0p**2+self.m0*sn_X_k0**2*sn_Y_k0**2
            DN_Z_k0 = DN_num / DN_den
            # step4 组合结果
            dsndZ = CN_Z_k0 * DN_Z_k0
            return dsndZ  
    
        # 初始化返回数组
        magV = np.zeros_like(Z, dtype=np.float64)
        u_vel = np.zeros_like(Z, dtype=np.float64)
        v_vel = np.zeros_like(Z, dtype=np.float64)

        zeta1 = SCT_g(np.real(self.z_so), np.imag(self.z_so))  # 源点变换
        zeta2 = SCT_g(np.real(self.z_si), np.imag(self.z_si))  # 汇点变换
    
        # 遍历 Z 数组中的每一个复数坐标
        it = np.nditer(Z, flags=['multi_index'])
        while not it.finished:
            z = it[0]
            zeta = SCT_g(np.real(z), np.imag(z))  # 当前点的 SC 变换
    
            # 计算 g'(Z)
            dsndz = SCT_g_prime(np.real(z), np.imag(z))

            # 打印调试信息
            #print(f"zeta:{zeta}, dsndz = {dsndz}")
    
            # 避免分母为零，添加 epsilon
            denominator_z1 = zeta - zeta1 + self.epsilon
            denominator_z2 = zeta - zeta2 + self.epsilon
    
            # 计算 dw/dZ
            dwdZ= (self.Q / np.pi) * dsndz * ((1/denominator_z1)-(1/denominator_z2))
    
            # 计算速度模、水平速度分量 u 和垂直速度分量 v
            magV[it.multi_index] = np.abs(dwdZ)  # 速度模
            u_vel[it.multi_index] = dwdZ.real  # 水平速度分量
            v_vel[it.multi_index] = -dwdZ.imag  # 垂直速度分量
    
            it.iternext()  # 移动到下一个坐标点
    
        return magV, u_vel, v_vel

    # 源汇压力
    def pressure(self, U_ref, magV, u_vel, v_vel):
        """
        计算压力分布和压力系数 CP。
    
        Parameters:
        ----------
        U_ref : float
            参考速度。
        magV : numpy.ndarray
            速度模，由 cplx_velocity 计算得到。
        u_vel : numpy.ndarray
            水平速度分量，由 cplx_velocity 计算得到。
        v_vel : numpy.ndarray
            垂直速度分量，由 cplx_velocity 计算得到。
    
        Returns:
        -------
        p : numpy.ndarray
            压力分布。
        CP : numpy.ndarray
            压力系数分布。
        d_v2 : numpy.ndarray
            速度误差项（用于诊断或调试）。
        """
        rho = 1.2  # 空气密度 (kg/m^3)
    
        # 计算速度平方的误差项 d_v2（用于诊断或调试）
        V_V = u_vel**2 + v_vel**2
        d_v2 = magV**2 - V_V
    
        # 计算压力分布 p 和压力系数 CP
        p = 0.5 * rho * (U_ref**2 - V_V)
        CP = p/(0.5*rho*U_ref**2)-0.99999
        #CP = 0.99999 - magV**2 / U_ref**2
    
        return p, CP, d_v2
    
    
# 室内气流绘图类: 绘制流函数、速度、压力分布
class RoomFlowPlot:
    def __init__(self, source_sink, Grid):
        """
        初始化 RoomFlowPlot 类
        
        Parameters:
        ----------
        X_size : list, tuple, or numpy array
            包含房间宽度范围的两个元素 [xmin, xmax]。
        Y_size : list, tuple, or numpy array
            包含房间高度范围的两个元素 [ymin, ymax]。
        source_sink : RoomSourceSink 对象
            包含源和汇的 RoomSourceSink 对象。
        Grid : Grid 对象
            网格对象，用于生成计算用的笛卡尔网格和复数网格。
        """
        
        # 对象属性化
        self.Grid = Grid  # 传入的网格对象
        self.plt_X, self.plt_Y = self.Grid.X, self.Grid.Y  # 笛卡尔网格坐标
        self.Z = self.Grid.to_complex()  # 复数网格坐标

        # 从 Grid 对象中获取 X 和 Y 范围
        self.plt_Xsize = 1.1 * self.Grid.X_ary[1]  # X 方向上比实际网格范围多 10%
        self.plt_Ysize = 1.1 * self.Grid.Y_ary[1]  # Y 方向上比实际网格范围多 10%
        
        # RoomSourceSink 对象属性化
        self.SoSi = source_sink  # 传入的 RoomSourceSink 对象

        # 计算 RoomSourceSink 的源点和汇点的复数坐标
        self.SoSi.compute_z_coordinates()

        # 设置源点和汇点的坐标
        self.plt_posX_source = np.real(self.SoSi.z_so)  # 源点的 x 坐标
        self.plt_posY_source = np.imag(self.SoSi.z_so)  # 源点的 y 坐标
        self.plt_posX_sink = np.real(self.SoSi.z_si)    # 汇点的 x 坐标
        self.plt_posY_sink = np.imag(self.SoSi.z_si)    # 汇点的 y 坐标

    def plot_real_geo(self, show_labels=False):
        """
        绘制z平面图像的源汇、壁面。
    
        Parameters:
        ----------
        show_labels : bool, optional
            是否显示坐标轴标签。默认显示。
        """
        #plt.figure(figsize=(8, 6))
        
        # 打印源和汇位置，确保它们在正确范围内
        print(f"Source Position: ({self.plt_posX_source}, {self.plt_posY_source})")
        print(f"Sink Position: ({self.plt_posX_sink}, {self.plt_posY_sink})")
    
        # 绘制源汇位置
        plt.scatter(self.plt_posX_source, self.plt_posY_source, color='r', s=80, marker="o", label='Source')
        plt.scatter(self.plt_posX_sink, self.plt_posY_sink, color='g', s=80, marker="X", label='Sink')
    
        # 绘制房间的墙壁（使用有限长度的线段）
        wall_left = self.plt_X[0][0]
        wall_right = self.plt_X[0][-1]
        wall_top = self.plt_Y[-1][0]
        wall_bottom = 0.0
    
        # 绘制顶部和底部墙壁（水平线段）
        plt.plot([wall_left, wall_right], [wall_top, wall_top], color='k', linestyle='-', linewidth=1)  # top wall
        plt.plot([wall_left, wall_right], [wall_bottom, wall_bottom], color='k', linestyle='-', linewidth=1)  # bottom wall
    
        # 绘制左侧和右侧墙壁（垂直线段）
        plt.plot([wall_left, wall_left], [wall_bottom, wall_top], color='k', linestyle='-', linewidth=1)  # left wall
        plt.plot([wall_right, wall_right], [wall_bottom, wall_top], color='k', linestyle='-', linewidth=1)  # right wall
        '''
        # 坐标轴范围
        plt.xlim([wall_left * 1.1, wall_right * 1.1])
        plt.ylim([wall_bottom, wall_top * 1.1])
    
        # 可选的标签和标题
        if show_labels:
            plt.xlabel('X Position (m)')
            plt.ylabel('Y Position (m)')
        plt.title('Room Geometry with Source and Sink')
    
        plt.legend()
        plt.grid(True)
        plt.show()
        '''
    def extract_contour_coordinates(self, phi_value):
        """
        从等值线绘图中提取给定等值线水平的坐标点。
    
        参数：
        contour_set: plt.contour 返回的 QuadContourSet 对象
        phi_value: 要提取的等值线水平值
    
        返回：
        coords: 包含该等值线所有坐标点的列表
        """
        coords = []
    
        # 直接遍历等值线对象的 paths 属性
        for path in self.contour_set.get_paths():
            # 获取该等值线的所有 (x, y) 坐标
            vertices = path.vertices
            # 对每个点进行检查，寻找与指定 phi 值最接近的坐标点
            if np.isclose(vertices[:, 0], phi_value).any() or np.isclose(vertices[:, 1], phi_value).any():
                coords.append(vertices)
    
        return coords

    
    def plot_psi(self, psi, title='Stream Function', num_levels=12):
        """
        绘制流函数等值线图。
    
        Parameters:
        ----------
        psi : numpy.ndarray
            流函数的二维数组，与网格对应。
        title : str, optional
            图形标题，默认是 'Stream Function'。
        num_levels : int, optional
            等值线的层数，默认是 12。
        """
        # Error handling for non-numeric values in psi
        if not np.issubdtype(psi.dtype, np.number):
            raise ValueError("psi must contain numeric values")
    
        # Calculate levels, ensuring they are valid
        min_psi = np.min(psi)
        max_psi = np.max(psi)
        if min_psi == max_psi:
            raise ValueError("psi must have a range of values for contouring")
        lvls = np.linspace(min_psi, max_psi, num_levels)

        # 动态设置图像大小，保持网格宽高比
        room_width = self.Grid.X_ary[1] - self.Grid.X_ary[0]
        room_height = self.Grid.Y_ary[1] - self.Grid.Y_ary[0]
        aspect_ratio = room_width / room_height

        # 设置figsize，保持合适比例
        plt.figure(figsize=(8, 10 / aspect_ratio))  # 根据网格宽高比动态调整图像大小
    
        # Create the plot
        #plt.figure(figsize=(8, 6))
        
        # 绘制几何图形，包括源汇位置和房间边界
        self.plot_real_geo()
    
        # Draw contours for psi (stream function)
        contour_psi = plt.contour(self.plt_X, self.plt_Y, psi, levels=lvls, cmap="coolwarm", linewidths=2)
        self.contour_set = contour_psi
        plt.clabel(contour_psi, inline=True, fontsize=10, fmt="%.2e")  # 使用科学计数法格式化标签

    
        # Add colorbar for better interpretation
        plt.colorbar(contour_psi, label='Stream Function')
    
        # Set plot limits to match the grid size
        plt.xlim(self.Grid.X_ary[0], self.Grid.X_ary[1])
        plt.ylim(self.Grid.Y_ary[0], self.Grid.Y_ary[1])
    
        # Labels and title
        plt.xlabel('X ')
        plt.ylabel('Y ')
        plt.title(title)
    
        # Grid settings
        plt.grid(True, linestyle='--', alpha=0.7)
    
        # Show the plot
        plt.show()


    def plot_phi(self, phi, title='Potential Function'):
        """
        绘制势函数的等值线图。
    
        Parameters:
        ----------
        phi : numpy.ndarray
            势函数的二维数组，与网格对应。
        title : str, optional
            图形标题，默认是 'Potential Function'。
        """
        # Error handling for non-numeric values in phi
        if not np.issubdtype(phi.dtype, np.number):
            raise ValueError("phi must contain numeric values")
    
        # 动态设置图像大小，保持网格宽高比
        room_width = self.Grid.X_ary[1] - self.Grid.X_ary[0]
        room_height = self.Grid.Y_ary[1] - self.Grid.Y_ary[0]
        aspect_ratio = room_width / room_height
    
        # 设置figsize，保持合适比例
        plt.figure(figsize=(8, 10 / aspect_ratio))  # 根据网格宽高比动态调整图像大小
        
        # Calculate levels for the contour plot
        min_phi = np.min(phi)
        max_phi = np.max(phi)
        if min_phi == max_phi:
            raise ValueError("phi must have a range of values for contouring")
        lvls = np.linspace(min_phi, max_phi, 12)
        
        # 势函数 phi 的等值线绘制
        contour_phi = plt.contour(self.plt_X, self.plt_Y, phi, levels=lvls, cmap="viridis", linewidths=2)
        plt.clabel(contour_phi, inline=True, fontsize=10, fmt="%.2e")  # 使用科学计数法格式化标签
        
        # 绘制几何图形，包括源汇位置和房间边界
        self.plot_real_geo()
        
        # 添加 colorbar 以便更好地解释势函数
        plt.colorbar(contour_phi, label='Potential Function')
    
        # Set plot limits to match the grid size
        plt.xlim(self.Grid.X_ary[0], self.Grid.X_ary[1])
        plt.ylim(self.Grid.Y_ary[0], self.Grid.Y_ary[1])
    
        # 设置坐标轴标签和标题
        plt.xlabel('X ')
        plt.ylabel('Y ')
        plt.title(title)
    
        # Grid settings
        plt.grid(True, linestyle='--', alpha=0.7)
    
        # Show the plot
        plt.show()


    def plot_flow_net(self, phi, psi, title='Flow Net'):
        """
        绘制流网图，即势函数 (phi) 和流函数 (psi) 的等值线图。
        
        Parameters:
        ----------
        phi : numpy.ndarray
            势函数的二维数组。
        psi : numpy.ndarray
            流函数的二维数组。
        title : str, optional
            图形标题，默认是 'Flow Net'。
        """
        # Error handling for non-numeric values in phi and psi
        if not np.issubdtype(phi.dtype, np.number) or not np.issubdtype(psi.dtype, np.number):
            raise ValueError("phi and psi must contain numeric values")
        
        # 动态设置图像大小，保持网格宽高比
        room_width = self.Grid.X_ary[1] - self.Grid.X_ary[0]
        room_height = self.Grid.Y_ary[1] - self.Grid.Y_ary[0]
        aspect_ratio = room_width / room_height
    
        # 设置figsize，保持合适比例
        plt.figure(figsize=(8, 10 / aspect_ratio))
        
        # Calculate levels for phi (Potential function) and psi (Stream function)
        min_phi, max_phi = np.min(phi), np.max(phi)
        min_psi, max_psi = np.min(psi), np.max(psi)
        
        if min_phi == max_phi or min_psi == max_psi:
            raise ValueError("phi and psi must have a range of values for contouring")
    
        lvls_phi = np.linspace(min_phi, max_phi, 12)
        lvls_psi = np.linspace(min_psi, max_psi, 12)
        
        # 绘制势函数 phi 的等值线
        contour_phi = plt.contour(self.plt_X, self.plt_Y, phi, levels=lvls_phi, cmap="viridis", linewidths=2)
        plt.clabel(contour_phi, inline=True, fontsize=10, fmt="%.2e")  # 使用科学计数法格式化标签

        
        # 绘制流函数 psi 的等值线
        contour_psi = plt.contour(self.plt_X, self.plt_Y, psi, levels=lvls_psi, cmap="coolwarm", linewidths=2)
        plt.clabel(contour_psi, inline=True, fontsize=10, fmt="%.2e")  # 使用科学计数法格式化标签

        
        # 调用散点绘制源汇位置和房间边界
        self.plot_real_geo()
    
        # Set plot limits to match the grid size
        plt.xlim(self.Grid.X_ary[0], self.Grid.X_ary[1])
        plt.ylim(self.Grid.Y_ary[0], self.Grid.Y_ary[1])
        
        # 添加坐标轴标签和标题
        plt.xlabel('X ')
        plt.ylabel('Y ')
        plt.title(title)
        plt.grid(True, linestyle='--', alpha=0.7)
        
        # 显示图像
        plt.show()


    def plot_vxvy(self, vx, vy, title='Velocity Field', density=1, linewidth=1):
        # 错误处理：确保 vx 和 vy 都是数值类型
        if not (np.issubdtype(vx.dtype, np.number) and np.issubdtype(vy.dtype, np.number)):
            raise ValueError("vx and vy must contain numeric values")
        
        # 确保数据有合理的范围
        if np.min(vx) == np.max(vx) or np.min(vy) == np.max(vy):
            raise ValueError("vx and vy must have a range of values for stream plotting")

        # 动态设置图像大小，保持网格宽高比
        room_width = self.Grid.X_ary[1] - self.Grid.X_ary[0]
        room_height = self.Grid.Y_ary[1] - self.Grid.Y_ary[0]
        aspect_ratio = room_width / room_height
    
        # 设置figsize，保持合适比例
        plt.figure(figsize=(8, 8 / aspect_ratio))
        
        # 绘制速度流线图
        plt.streamplot(self.plt_X, self.plt_Y, vx, vy, density=density, linewidth=linewidth, arrowsize=1, arrowstyle='->')
        
        # 调用函数绘制源汇位置
        self.plot_real_geo()  
        
        # 图像美化：添加标签、网格等
        plt.xlabel('x')
        plt.ylabel('y')
        plt.title(title)
        plt.grid(True, linestyle='--', alpha=0.7)
        
        # 显示图像
        plt.show()


    def plot_MagV(self, magV, title='Velocity Magnitude'):
        # 动态设置图像大小，保持网格宽高比
        room_width = self.Grid.X_ary[1] - self.Grid.X_ary[0]
        room_height = self.Grid.Y_ary[1] - self.Grid.Y_ary[0]
        aspect_ratio = room_width / room_height
    
        # 设置figsize，保持合适比例
        plt.figure(figsize=(8, 10 / aspect_ratio))

        contour_magV = plt.contour(self.plt_X, self.plt_Y, magV, cmap="viridis", linewidths=2)   
        #plt.clabel(contour_magV, inline=True, fontsize=10)
        self.plot_real_geo()  # 调用散点绘制源汇位置
        plt.colorbar(contour_magV, label='Velocity Magnitude')
        plt.xlabel('Width')
        plt.ylabel('Height')
        plt.title(title)
        plt.grid(True)
        plt.show()

    def plot_pressure(self, p, title='Pressure Coefficient Distribution'):
        # 动态设置图像大小，保持网格宽高比
        room_width = self.Grid.X_ary[1] - self.Grid.X_ary[0]
        room_height = self.Grid.Y_ary[1] - self.Grid.Y_ary[0]
        aspect_ratio = room_width / room_height
    
        # 设置figsize，保持合适比例
        plt.figure(figsize=(8, 10 / aspect_ratio))
        #lvlsp = np.array([np.min(self.plt_p), np.max(self.plt_p)])
        #print(f"p:{lvlsp}")
        contour = plt.contour(self.plt_X, self.plt_Y, p,
                             levels=20, cmap=plt.cm.viridis, linewidths=1)
        #plt.contourf(self.room_flow.room.X, self.room_flow.room.Y, self.plt_p, cmap=plt.cm.plasma)
        plt.clabel(contour, inline=True, fontsize=10)
        #self.plot_real_geo()  # 调用散点绘制源汇位置
        plt.colorbar(contour, label='Pressure Coefficient(Cp)',extend="both")
        plt.xlabel('x')
        plt.ylabel('y')
        plt.title(title)
        plt.grid(True)
        plt.show()
    
# 调用
def calculate_and_plot_geo_domain(length, height, v_inlet, d, L1, H1, L2, H2, N=50):
    """
    计算和绘制几何坐标系（无量纲）下的流场特性。

    Parameters:
    -----------
    length : float
        房间长度（物理量）。
    height : float
        房间高度（物理量）。
    v_inlet : float
        入口速度。
    d : float
        源汇强度的尺度因子。
    L1, H1 : float
        源的位置 (x, y) 坐标。
    L2, H2 : float
        汇的位置 (x, y) 坐标。
    N : int, optional
        网格划分数量，默认值为 50。
    """
    # 几何域（无量纲）参数
    K_p = height / height
    K2 = length / height
    geo_x = np.array([-K2 / 2, K2 / 2])
    geo_y = np.array([0.0, K_p])
    x1 = length
    y1 = height
    room_geo = Grid(length, height, N)

    # SC变换参数计算：C0, k0
    SC_solver = ParameterSolver(length, height)
    C0, k0, _ = SC_solver.solve()

    # 源汇强度
    Strength_sosi = v_inlet * d
    room_sosi_geo = RoomSourceSink(length, height, L1, H1, L2, H2, Strength_sosi, C0, k0)

    # 计算复位势、速度场和压力分布
    W_sosi_geo, phi_sosi_geo, psi_sosi_geo = room_sosi_geo.cplx_W(room_geo.X, room_geo.Y)
    V_sosi_geo, vx_sosi_geo, vy_sosi_geo = room_sosi_geo.cplx_velocity(room_geo.X, room_geo.Y)
    p_sosi_geo, CP_sosi_geo, _ = room_sosi_geo.pressure(v_inlet, V_sosi_geo, vx_sosi_geo, vy_sosi_geo)

    # 绘图
    plt_sosi_geo = RoomFlowPlot(room_sosi_geo, room_geo)
    plt_sosi_geo.plot_psi(psi_sosi_geo, title='Stream Function (Dimensionless)')
    plt_sosi_geo.plot_phi(phi_sosi_geo, title='Potential Function (Dimensionless)')
    plt_sosi_geo.plot_flow_net(phi_sosi_geo, psi_sosi_geo, title='Flow Net (Dimensionless)')
    plt_sosi_geo.plot_vxvy(vx_sosi_geo, vy_sosi_geo, title='Velocity Field (Dimensionless)')
    plt_sosi_geo.plot_MagV(V_sosi_geo, title='Velocity Magnitude (Dimensionless)')
    plt_sosi_geo.plot_pressure(CP_sosi_geo)


def calculate_and_plot_phy_domain(length, height, v_inlet, d, L1, H1, L2, H2, N=5):
    """
    计算和绘制物理坐标系（有量纲）下的流场特性。

    Parameters:
    -----------
    length : float
        房间长度（物理量）。
    height : float
        房间高度（物理量）。
    v_inlet : float
        入口速度。
    d : float
        源汇强度的尺度因子。
    L1, H1 : float
        源的位置 (x, y) 坐标。
    L2, H2 : float
        汇的位置 (x, y) 坐标。
    N : int, optional
        网格划分数量，默认值为 5。
    """
    # 物理域参数
    phy_x = np.array([-length / 2, length / 2])
    phy_y = np.array([0.0, height])
    room_phy = Grid(phy_x, phy_y, N)

    # SC变换参数计算：C0, k0
    SC_solver = ParameterSolver(length, height)
    C0, k0, _ = SC_solver.solve()

    # 源汇强度
    Strength_sosi = v_inlet * d
    room_sosi_phy = RoomSourceSink(length, height, L1, H1, L2, H2, Strength_sosi, C0, k0)

    # 计算复位势、速度场和压力分布
    W_sosi_phy, phi_sosi_phy, psi_sosi_phy = room_sosi_phy.cplx_W(room_phy.X, room_phy.Y)
    V_sosi_phy, vx_sosi_phy, vy_sosi_phy = room_sosi_phy.cplx_velocity(room_phy.X, room_phy.Y)
    p_sosi_phy, CP_sosi_phy, _ = room_sosi_phy.pressure(v_inlet, V_sosi_phy, vx_sosi_phy, vy_sosi_phy)

    # 绘图
    plt_sosi_phy = RoomFlowPlot(room_sosi_phy, room_phy)
    plt_sosi_phy.plot_psi(psi_sosi_phy, title='Stream Function (Physical Domain)')
    plt_sosi_phy.plot_phi(phi_sosi_phy, title='Potential Function (Physical Domain)')
    plt_sosi_phy.plot_flow_net(phi_sosi_phy, psi_sosi_phy, title='Flow Net (Physical Domain)')
    plt_sosi_phy.plot_vxvy(vx_sosi_phy, vy_sosi_phy, title='Velocity Field (Physical Domain)')
    plt_sosi_phy.plot_MagV(V_sosi_phy, title='Velocity Magnitude (Physical Domain)')
    plt_sosi_phy.plot_pressure(CP_sosi_phy, title='Pressure Field (Physical Domain)')


class SourceSink2:
    """
    定义室内一源一汇流动（Source-sink Flow）的基本计算。

    该类用于计算势流的复势函数、流函数、速度分布等参数。
    """
    def __init__(self, sosi_strength, source_org, sink_org, C, k, length, height):
        """
        初始化一源一汇流动参数

        Parameters
        ----------
        sosi_strength : float
            Strength of the source sink pair (m^2/s)
        source_org : tuple or complex
            Source origin coordinates (length_1, height_1) in meters
        sink_org : tuple or complex
            Sink origin coordinates (length_2, height_2) in meters
        C : float
            Parameter for the Schwarz-Christoffel (SC) transformation function
        k : float
            Modulus number for elliptic functions
        length : float
            Length of room (m)
        height : float
            Height of room (m)
        """
        # 源和汇的强度
        self.Q_source = sosi_strength[0]
        self.Q_sink = sosi_strength[1]

        # 房间尺寸
        self.L = length
        self.H = height

        # 源和汇的坐标
        if isinstance(source_org, (tuple, list)) and len(source_org) == 2:
            self.L1, self.H1 = source_org
        elif isinstance(source_org, complex):
            self.L1, self.H1 = source_org.real, source_org.imag
        else:
            raise ValueError("source_org must be a tuple/list (x, y) or complex number")
        if isinstance(sink_org, (tuple, list)) and len(sink_org) == 2:
            self.L2, self.H2 = sink_org
        elif isinstance(sink_org, complex):
            self.L2, self.H2 = sink_org.real, sink_org.imag
        else:
            raise ValueError("sink_org must be a tuple/list (x, y) or complex number")

        # 检查坐标是否在房间范围内（移到 __init__ 中统一验证）
        if not (0 <= self.L1 <= self.L and 0 <= self.H1 <= self.H):
            raise ValueError("Source coordinates must be within room dimensions")
        if not (0 <= self.L2 <= self.L and 0 <= self.H2 <= self.H):
            raise ValueError("Sink coordinates must be within room dimensions")

        # 初始化无量纲坐标
        self.z_so = None
        self.z_si = None

        # SC 变换参数
        self.C0 = C
        self.k0 = k
        self.m0 = k ** 2
        self.k0p = np.sqrt(1 - k ** 2)
        self.m0p = 1 - self.m0
        self.epsilon = 1e-10

    def __repr__(self):
        """
        返回类的字符串表示，方便打印查看
        """
        return (f"SourceSink(源强={self.Q_source}, 汇强={self.Q_sink}, "
                f"源坐标=({self.L1}, {self.H1}), 汇坐标=({self.L2}, {self.H2}), "
                f"C={self.C0}, k={self.k0}, length={self.L}, height={self.H})")

    def compute_z_coordinates(self):
        '''
        计算源点和汇点的无量纲坐标，并保存到类变量 z_so 和 z_si 中
        '''
        E = self.L1 - self.L / 2 + 1j * self.H1  # 源点
        F = self.L2 - self.L / 2 + 1j * self.H2  # 汇点
        self.z_so = E / self.H  # 源点无量纲坐标
        self.z_si = F / self.H  # 汇点无量纲坐标

        # print(f"计算得到的源点坐标 z_so: {self.z_so}")
        # print(f"计算得到的汇点坐标 z_si: {self.z_si}")

    # 源汇复势
    def cplx_W(self, X, Y):
        """
        计算复位势及其分量，输入X,Y坐标均为无量纲坐标
        Parameters:
            X : np.ndarray
                Cartesian X coordinates
            Y : np.ndarray
                Cartesian Y coordinates
        Returns:
            w_ptn : np.ndarray
                Complex potential
            phi : np.ndarray
                Potential function
            psi : np.ndarray
                Stream function
        """
        # 确保先计算源点和汇点的复数坐标
        if self.z_so is None or self.z_si is None:
            self.compute_z_coordinates()

        Z = X + Y * 1j  # 复数形式的坐标 Z = X + iY

        # 定义SC变换函数 zeta = sn(z / C0; k**2)
        def SCT_g(X, Y):
            """
            计算复数域的椭圆正弦函数 sn(z, k)
            Parameters:
            X -- 实部
            Y -- 虚部
            Returns:
            复数域的椭圆正弦函数结果 zeta
            """
            # 归一化
            use_X = X / self.C0
            use_Y = Y / self.C0

            # 椭圆函数的计算
            sn_X_k0, cn_X_k0, dn_X_k0, _ = sp.ellipj(use_X, self.m0)
            sn_Y_k0, cn_Y_k0, dn_Y_k0, _ = sp.ellipj(use_Y, self.m0)

            sn_X_k0p, cn_X_k0p, dn_X_k0p, _ = sp.ellipj(use_X, self.m0p)
            sn_Y_k0p, cn_Y_k0p, dn_Y_k0p, _ = sp.ellipj(use_Y, self.m0p)

            # 计算分子和分母
            num = sn_X_k0 * dn_Y_k0p + 1j * sn_Y_k0p * cn_Y_k0p * cn_X_k0 * dn_X_k0
            den = cn_Y_k0p ** 2 + self.k0 ** 2 * sn_X_k0 ** 2 * sn_Y_k0p ** 2

            zeta = num / den
            return zeta

        # 初始化返回数组
        w_ptn = np.zeros_like(Z, dtype=np.complex128)
        phi = np.zeros_like(Z, dtype=np.float64)
        psi = np.zeros_like(Z, dtype=np.float64)

        zeta1 = SCT_g(np.real(self.z_so), np.imag(self.z_so))  # 源点变换
        zeta2 = SCT_g(np.real(self.z_si), np.imag(self.z_si))  # 汇点变换

        # 遍历 Z 数组
        it = np.nditer(Z, flags=['multi_index'])
        while not it.finished:
            z = it[0]
            zeta = SCT_g(np.real(z), np.imag(z))  # 当前点的 SC 变换

            # 避免分母为零，添加 epsilon
            # denominator_z1 = zeta - zeta1 + self.epsilon
            # denominator_z2 = zeta - zeta2 + self.epsilon

            # 计算复势函数 w_ptn
            w_ptn[it.multi_index] = (self.Q / np.pi) * (np.log(zeta - zeta1) - np.log(zeta - zeta2))

            # 提取位势和流函数
            phi[it.multi_index] = w_ptn[it.multi_index].real
            psi[it.multi_index] = w_ptn[it.multi_index].imag

            it.iternext()  # 移动到下一个坐标点

        return w_ptn, phi, psi

    # 源汇复速度
    def cplx_velocity(self, X, Y):
        """
        计算复势 w 对复坐标 Z 的导数
        :param X: 笛卡尔坐标 X
        :param Y: 笛卡尔坐标 Y
        :return: dw_dZ (复势的导数)
        """
        # 确保先计算源点和汇点的复数坐标
        if self.z_so is None or self.z_si is None:
            self.compute_z_coordinates()

        Z = X + Y * 1j  # 复数形式的坐标 Z = X + iY

        # 定义SC变换函数 zeta = sn(z / C0; k**2)
        def SCT_g(X, Y):
            """
            计算复数域的椭圆正弦函数 sn(z, k)
            Parameters:
            X -- 实部
            Y -- 虚部
            Returns:
            复数域的椭圆正弦函数结果 zeta
            """
            # 归一化
            use_X = X / self.C0
            use_Y = Y / self.C0

            # 椭圆函数的计算
            sn_X_k0, cn_X_k0, dn_X_k0, _ = sp.ellipj(use_X, self.m0)
            sn_Y_k0, cn_Y_k0, dn_Y_k0, _ = sp.ellipj(use_Y, self.m0)

            sn_X_k0p, cn_X_k0p, dn_X_k0p, _ = sp.ellipj(use_X, self.m0p)
            sn_Y_k0p, cn_Y_k0p, dn_Y_k0p, _ = sp.ellipj(use_Y, self.m0p)

            # 计算分子和分母
            num = sn_X_k0 * dn_Y_k0p + 1j * sn_Y_k0p * cn_Y_k0p * cn_X_k0 * dn_X_k0
            den = cn_Y_k0p ** 2 + self.k0 ** 2 * sn_X_k0 ** 2 * sn_Y_k0p ** 2

            zeta = num / den
            return zeta

        # 椭圆正弦函数的导数 dSCTdZ(z / C0; k**2)
        def SCT_g_prime(X, Y):
            """
            z = x + iy
            计算椭圆正弦函数的导数 cn(z, k) * dn(z, k)
            Z -- 复数形式的坐标
            """
            # step0 归一化
            use_X = X / self.C0
            use_Y = Y / self.C0
            # step1 实部、虚部分离计算模数的椭圆函数
            sn_X_k0, cn_X_k0, dn_X_k0, _ = sp.ellipj(use_X, self.m0)
            sn_Y_k0, cn_Y_k0, dn_Y_k0, _ = sp.ellipj(use_Y, self.m0)
            # step2 实部、虚部分离计算补模数的椭圆函数
            sn_X_k0p, cn_X_k0p, dn_X_k0p, _ = sp.ellipj(use_X, self.m0p)
            sn_Y_k0p, cn_Y_k0p, dn_Y_k0p, _ = sp.ellipj(use_Y, self.m0p)
            # step3 组件计算
            sn_iY_k0 = 1j * (sn_Y_k0p / cn_Y_k0p)
            cn_iY_k0 = 1 / cn_Y_k0p
            dn_iY_k0 = dn_Y_k0p / cn_Y_k0p

            CN_num = sn_X_k0 * cn_X_k0 * dn_iY_k0 - sn_iY_k0 * cn_iY_k0 * dn_X_k0
            CN_den = sn_X_k0 * cn_iY_k0 * dn_iY_k0 - sn_iY_k0 * cn_X_k0 * dn_X_k0
            CN_Z_k0 = CN_num / CN_den

            DN_num = dn_X_k0 * cn_Y_k0p * dn_Y_k0p - 1j * self.m0 * sn_X_k0 * cn_X_k0 * sn_Y_k0p
            DN_den = cn_Y_k0p ** 2 + self.m0 * sn_X_k0 ** 2 * sn_Y_k0 ** 2
            DN_Z_k0 = DN_num / DN_den
            # step4 组合结果
            dsndZ = CN_Z_k0 * DN_Z_k0
            return dsndZ

            # 初始化返回数组

        magV = np.zeros_like(Z, dtype=np.float64)
        u_vel = np.zeros_like(Z, dtype=np.float64)
        v_vel = np.zeros_like(Z, dtype=np.float64)

        zeta1 = SCT_g(np.real(self.z_so), np.imag(self.z_so))  # 源点变换
        zeta2 = SCT_g(np.real(self.z_si), np.imag(self.z_si))  # 汇点变换

        # 遍历 Z 数组中的每一个复数坐标
        it = np.nditer(Z, flags=['multi_index'])
        while not it.finished:
            z = it[0]
            zeta = SCT_g(np.real(z), np.imag(z))  # 当前点的 SC 变换

            # 计算 g'(Z)
            dsndz = SCT_g_prime(np.real(z), np.imag(z))

            # 打印调试信息
            # print(f"zeta:{zeta}, dsndz = {dsndz}")

            # 避免分母为零，添加 epsilon
            denominator_z1 = zeta - zeta1 + self.epsilon
            denominator_z2 = zeta - zeta2 + self.epsilon

            # 计算 dw/dZ
            dwdZ = (self.Q / np.pi) * dsndz * ((1 / denominator_z1) - (1 / denominator_z2))

            # 计算速度模、水平速度分量 u 和垂直速度分量 v
            magV[it.multi_index] = np.abs(dwdZ)  # 速度模
            u_vel[it.multi_index] = dwdZ.real  # 水平速度分量
            v_vel[it.multi_index] = -dwdZ.imag  # 垂直速度分量

            it.iternext()  # 移动到下一个坐标点

        return magV, u_vel, v_vel

    # 源汇压力
    def pressure(self, U_ref, magV, u_vel, v_vel):
        """
        计算压力分布和压力系数 CP。

        Parameters:
        ----------
        U_ref : float
            参考速度。
        magV : numpy.ndarray
            速度模，由 cplx_velocity 计算得到。
        u_vel : numpy.ndarray
            水平速度分量，由 cplx_velocity 计算得到。
        v_vel : numpy.ndarray
            垂直速度分量，由 cplx_velocity 计算得到。

        Returns:
        -------
        p : numpy.ndarray
            压力分布。
        CP : numpy.ndarray
            压力系数分布。
        d_v2 : numpy.ndarray
            速度误差项（用于诊断或调试）。
        """
        rho = 1.2  # 空气密度 (kg/m^3)

        # 计算速度平方的误差项 d_v2（用于诊断或调试）
        V_V = u_vel ** 2 + v_vel ** 2
        d_v2 = magV ** 2 - V_V

        # 计算压力分布 p 和压力系数 CP
        p = 0.5 * rho * (U_ref ** 2 - V_V)
        CP = p / (0.5 * rho * U_ref ** 2) - 0.99999
        # CP = 0.99999 - magV**2 / U_ref**2

        return p, CP, d_v2

"""____示例____"""
# Example of using the general models
if __name__ == "__main__":

    '''
    # 定义物理参数
    length = 5.0  # 房间长度 (m)
    height = 3.0  # 房间高度 (m)

    v_inlet = 40  # 入口速度 (m/s)

    d = 0.05      # 源汇强度的尺度因子

    L1, H1 = 0.1, height  # 源位置(m)
    L2, H2 = length - 0.0, 1.4  # 汇位置(m)

    # 使用几何域计算和绘图
    calculate_and_plot_geo_domain(length, height, v_inlet, d, L1, H1, L2, H2)
    
    # (TODO)使用物理域计算和绘图
    #calculate_and_plot_phy_domain(length, height, v_inlet, d, L1, H1, L2, H2)
    '''
    # 计算域物理参数(有量纲)
    height = 3.0 # m
    length = 5.0 # m
    phy_x = np.array([-length/2, length/2])
    phy_y = np.array([0.0, height])
    N = 5
    room_phy = Grid(length, height, N) # 网格坐标
    # 计算域几何参数(无量纲)
    K_p = height / height
    K2 = length / height
    geo_x = np.array([-K2/2, K2/2])
    geo_y = np.array([0.0, K_p])
    N = 50
    room_geo = Grid(K_p, K2, N) # 网格坐标
    
    # SC变换参数计算：C0,k0
    SC_solver = ParameterSolver(length, height)
    C0, k0, k0_prime = SC_solver.solve()# 计算C0模数 k0 和补模数 k0'
    print(f"计算得到的模数 k0 = {k0}")
    print(f"补模数 k0' = {k0_prime}")
    print(f"计算得到的系数 C0 = {C0}")
    
    # 源汇的物理参数
    v_inlet = 40 # m/s
    d = 0.05 # m
    Strength_sosi = v_inlet * d
    # 源汇物理位置距离(m)
    L1 = 0.1
    H1 = height
    L2 = length - 0.1
    H2 = height
    
    # 源汇模型
    room_sosi = RoomSourceSink(length, height, L1, H1, L2, H2, Strength_sosi, C0, k0)
    W_sosi, phi_sosi, psi_sosi = room_sosi.cplx_W(room_geo.X, room_geo.Y)# 复位势
    #print(f"phi:{phi_sosi}")
    #print(f"psi:{psi_sosi}")
    # 复速度
    V_sosi, vx_sosi, vy_sosi = room_sosi.cplx_velocity(room_geo.X, room_geo.Y)
    # 压力分布
    p_sosi, CP_sosi,dv2 = room_sosi.pressure(v_inlet, V_sosi, vx_sosi, vy_sosi)
    #print(dv2)
    #print(CP_sosi)
    #print(vy_sosi)
    
    # 绘图
    #target_psi = 0.0182
    plt_sosi = RoomFlowPlot(room_sosi, room_geo)
    plt_sosi.plot_psi(psi_sosi, title='Stream Function')
    #print(psi_sosi)
    
    plt_sosi.plot_phi(phi_sosi, title='Potential Function')
    plt_sosi.plot_flow_net(phi_sosi, psi_sosi, title='Flow Net')
    plt_sosi.plot_vxvy(vx_sosi, vy_sosi, title='Velocity Field')
    
    #plt_sosi.plot_vxvy(RR, II, title='2 order')
    plt_sosi.plot_MagV(V_sosi, title='Velocity Magnitude')
    plt_sosi.plot_pressure(CP_sosi, title='Pressure Field')
